#include "CellParticles.h"
#include "../Input/PicParams.h"

#include <math.h>
#include <iostream>
#include <algorithm>


template<int dim>
CellParticles<dim>::CellParticles()
{

}

template<int dim>
void CellParticles<dim>::initialize(const PicParams* params, const SpeciesParameter& species_param, const vector<double>& cell_min, const vector<double>& cell_length)
{
    int n_particles_reserve = species_param.c_part_max * species_param.n_part_per_cell_for_weight;
    data_.reserve(n_particles_reserve);

    int n_particles = species_param.n_part_per_cell;
    data_.resize(n_particles);

    //std::cout<<"n_particles cellparticle init: "<<n_particles<<endl;

    for(int i_particle = 0; i_particle < n_particles; i_particle++)
    {
        //initialize velocity
        if (species_param.init_velocity_type == "cold")
        {
            for(int i = 0; i < 3; i++)
            {
                data_[i_particle].velocity[i] = 0.0 + species_param.mean_velocity[i];
            }
        }
        else if(species_param.init_velocity_type == "maxwell")
        {
            if(!params->is_temperature_anisotropic)
            {
                double vt = sqrt(2.0 * species_param.temperature * params->const_e / species_param.mass);
                double r1;
                double r2;
                for(int i = 0; i < 3; i++)
                {
                    do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

                    r2 = (double)rand() / RAND_MAX;
                    data_[i_particle].velocity[i] = vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2) + species_param.mean_velocity[i];
                }
            }
            else
            {
                double temperature[3];
                temperature[0] = species_param.temperature_x;
                temperature[1] = species_param.temperature_y;
                temperature[2] = species_param.temperature_z;
                for(int i = 0; i < 3; i++)
                {
                    double vt = sqrt(2.0 * temperature[i] * params->const_e / species_param.mass);
                    double r1;
                    double r2;
                    do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

                    r2 = (double)rand() / RAND_MAX;
                    data_[i_particle].velocity[i] = vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2) + species_param.mean_velocity[i];
                }
            }

        }

        //initialize position
        for(int i = 0; i < dim; i++) 
        {
            if(species_param.init_position_type == "regular")
            {
                data_[i_particle].position[i] = cell_min[i] + (i_particle + 0.5) * cell_length[i] / n_particles;
            } 
            else if(species_param.init_position_type == "random")
            {
                data_[i_particle].position[i] = cell_min[i] + (((double)rand() / RAND_MAX)) * cell_length[i];
            }
        }
    }    
    
}


template<int dim>
void CellParticles<dim>::reinitialize_particles(const PicParams* params, const SpeciesParameter& species_param, const vector<double>& cell_min,
                                                const vector<double>& cell_length, const vector<double>& mean_velocity, double temperature, int number_particles)
{
    int n_particles = number_particles;
    data_.resize(n_particles);

    //std::cout<<"n_particles cellparticle init: "<<n_particles<<endl;

    for(int i_particle = 0; i_particle < n_particles; i_particle++)
    {
        //initialize velocity
        if (species_param.init_velocity_type == "cold")
        {
            for(int i = 0; i < 3; i++)
            {
                data_[i_particle].velocity[i] = 0.0 + mean_velocity[i];
            }
        }
        else if(species_param.init_velocity_type == "maxwell")
        {
            double vt = sqrt(2.0 * temperature * params->const_e / species_param.mass);
            double r1;
            double r2;
            for(int i = 0; i < 3; i++)
            {
                do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

                r2 = (double)rand() / RAND_MAX;
                data_[i_particle].velocity[i] = vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2) + mean_velocity[i];
            }
        }

        //initialize position
        for(int i = 0; i < dim; i++)
        {
            if(species_param.init_position_type == "regular")
            {
                data_[i_particle].position[i] = cell_min[i] + (i_particle + 0.5) * cell_length[i] / n_particles;
            }
            else if(species_param.init_position_type == "random")
            {
                data_[i_particle].position[i] = cell_min[i] + (((double)rand() / RAND_MAX)) * cell_length[i];
            }
        }
    }

}

template<int dim>
void CellParticles<dim>::add_particles(const PicParams* params, const SpeciesParameter& species_param, const vector<double>& cell_min,
                                                const vector<double>& cell_length, const vector<double>& mean_velocity, double temperature, int number_particles)
{
    int n_particles_pre = data_.size();
    data_.resize(n_particles_pre + number_particles);

    //std::cout<<"n_particles cellparticle init: "<<n_particles<<endl;

    for(int i_particle = n_particles_pre; i_particle < data_.size(); i_particle++)
    {
        //initialize velocity
        if (species_param.init_velocity_type == "cold")
        {
            for(int i = 0; i < 3; i++)
            {
                data_[i_particle].velocity[i] = 0.0 + mean_velocity[i];
            }
        }
        else if(species_param.init_velocity_type == "maxwell")
        {
            double vt = sqrt(2.0 * temperature * params->const_e / species_param.mass);
            double r1;
            double r2;
            for(int i = 0; i < 3; i++)
            {
                do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

                r2 = (double)rand() / RAND_MAX;
                data_[i_particle].velocity[i] = vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2) + mean_velocity[i];
            }
        }

        //initialize position
        for(int i = 0; i < dim; i++)
        {
            if(species_param.init_position_type == "regular")
            {
                data_[i_particle].position[i] = cell_min[i] + (i_particle - n_particles_pre + 0.5) * cell_length[i] / number_particles;
            }
            else if(species_param.init_position_type == "random")
            {
                data_[i_particle].position[i] = cell_min[i] + (((double)rand() / RAND_MAX)) * cell_length[i];
            }
        }
    }

}



template<int dim>
void CellParticles<dim>::add_particles(int number_particles)
{
    int n_particles_pre = data_.size();
    data_.resize(n_particles_pre + number_particles);
}


template<int dim>
void CellParticles<dim>::heat(const PicParams* params, const SpeciesParameter& species_param, double temperature_heat)
{
    double velocity_unit[3];
    double velocity_magnitude;
    double velocity_magnitude_heat;

    for(int i_particle = 0; i_particle < data_.size(); i_particle++)
    {
        velocity_magnitude = sqrt( data_[i_particle].velocity[0] * data_[i_particle].velocity[0]
                                 + data_[i_particle].velocity[1] * data_[i_particle].velocity[1]
                                 + data_[i_particle].velocity[2] * data_[i_particle].velocity[2] );
        if(velocity_magnitude == 0.0)
        {
          continue;
        }
        velocity_unit[0] = data_[i_particle].velocity[0] / velocity_magnitude;
        velocity_unit[1] = data_[i_particle].velocity[1] / velocity_magnitude;
        velocity_unit[2] = data_[i_particle].velocity[2] / velocity_magnitude;

        velocity_magnitude_heat = sqrt(2.0 * temperature_heat * params->const_e / species_param.mass + velocity_magnitude * velocity_magnitude);

        data_[i_particle].velocity[0] = velocity_unit[0] * velocity_magnitude_heat;
        data_[i_particle].velocity[1] = velocity_unit[1] * velocity_magnitude_heat;
        data_[i_particle].velocity[2] = velocity_unit[2] * velocity_magnitude_heat;
    }

}


template<int dim>
void CellParticles<dim>::heat_new(const PicParams* params, const SpeciesParameter& species_param, const ParticleSourceParameter& ps_param, double energy_heat)
{
    double velocity_unit[3];
    double velocity_magnitude;
    double velocity_magnitude_heat;

    double number_particle_maxwell_double = energy_heat / (1.5 * ps_param.temperature);
    int number_particle_maxwell = number_particle_maxwell_double;
    double energy_rem = (number_particle_maxwell_double - number_particle_maxwell) * ps_param.temperature;

    if(number_particle_maxwell > data_.size())
    {
        log_fatal<<"Error in cell particles heat_new(): "<<number_particle_maxwell<<" "<<data_.size();
        return;
    }

    for(int i_particle = 0; i_particle < number_particle_maxwell; i_particle++)
    {

        double vt = sqrt(2.0 * ps_param.temperature * params->const_e / species_param.mass);
        double r1;
        double r2;
        double velocity_magnitude_add = 0.0;
        for(int i = 0; i < 3; i++)
        {
            do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

            r2 = (double)rand() / RAND_MAX;
            velocity_magnitude_add += vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2);
        }



        velocity_magnitude = sqrt( data_[i_particle].velocity[0] * data_[i_particle].velocity[0]
                                 + data_[i_particle].velocity[1] * data_[i_particle].velocity[1]
                                 + data_[i_particle].velocity[2] * data_[i_particle].velocity[2] );
        if(velocity_magnitude == 0.0)
        {
          continue;
        }
        velocity_unit[0] = data_[i_particle].velocity[0] / velocity_magnitude;
        velocity_unit[1] = data_[i_particle].velocity[1] / velocity_magnitude;
        velocity_unit[2] = data_[i_particle].velocity[2] / velocity_magnitude;

        velocity_magnitude_heat = sqrt(velocity_magnitude_add * velocity_magnitude_add + velocity_magnitude * velocity_magnitude);

        data_[i_particle].velocity[0] = velocity_unit[0] * velocity_magnitude_heat;
        data_[i_particle].velocity[1] = velocity_unit[1] * velocity_magnitude_heat;
        data_[i_particle].velocity[2] = velocity_unit[2] * velocity_magnitude_heat;
    }

    int i_final;
    if(number_particle_maxwell < data_.size())
    {
        i_final = number_particle_maxwell;
    }
    else
    {
        i_final = 0;
    }

    velocity_magnitude = sqrt( data_[i_final].velocity[0] * data_[i_final].velocity[0]
                             + data_[i_final].velocity[1] * data_[i_final].velocity[1]
                             + data_[i_final].velocity[2] * data_[i_final].velocity[2] );
    if(velocity_magnitude == 0.0)
    {
      return;
    }
    velocity_unit[0] = data_[i_final].velocity[0] / velocity_magnitude;
    velocity_unit[1] = data_[i_final].velocity[1] / velocity_magnitude;
    velocity_unit[2] = data_[i_final].velocity[2] / velocity_magnitude;

    velocity_magnitude_heat = sqrt(2.0 * energy_rem * params->const_e / species_param.mass + velocity_magnitude * velocity_magnitude);

    data_[i_final].velocity[0] = velocity_unit[0] * velocity_magnitude_heat;
    data_[i_final].velocity[1] = velocity_unit[1] * velocity_magnitude_heat;
    data_[i_final].velocity[2] = velocity_unit[2] * velocity_magnitude_heat;

}


template<int dim>
void CellParticles<dim>::heat_new2(const PicParams* params, const SpeciesParameter& species_param, const ParticleSourceParameter& ps_param, 
                                   double energy_heat, double heat_number_ratio)
{
    double velocity_unit[3];
    double velocity_magnitude;
    double velocity_magnitude_heat;

    double temperature_heat = 2.0 * energy_heat / 3.0;

    int heat_number = data_.size() * heat_number_ratio;
    if(heat_number == 0)
    {
        heat_number = 1;
    }

    temperature_heat /= heat_number;

    for(int i_particle = 0; i_particle < heat_number; i_particle++)
    {

        double vt = sqrt(2.0 * temperature_heat * params->const_e / species_param.mass);
        double r1;
        double r2;
        double velocity_magnitude_add = 0.0;
        for(int i = 0; i < 3; i++)
        {
            do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

            r2 = (double)rand() / RAND_MAX;
            velocity_magnitude_add += vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2);
        }



        velocity_magnitude = sqrt( data_[i_particle].velocity[0] * data_[i_particle].velocity[0]
                                 + data_[i_particle].velocity[1] * data_[i_particle].velocity[1]
                                 + data_[i_particle].velocity[2] * data_[i_particle].velocity[2] );
        if(velocity_magnitude == 0.0)
        {
          continue;
        }
        velocity_unit[0] = data_[i_particle].velocity[0] / velocity_magnitude;
        velocity_unit[1] = data_[i_particle].velocity[1] / velocity_magnitude;
        velocity_unit[2] = data_[i_particle].velocity[2] / velocity_magnitude;

        velocity_magnitude_heat = sqrt(velocity_magnitude_add * velocity_magnitude_add + velocity_magnitude * velocity_magnitude);

        data_[i_particle].velocity[0] = velocity_unit[0] * velocity_magnitude_heat;
        data_[i_particle].velocity[1] = velocity_unit[1] * velocity_magnitude_heat;
        data_[i_particle].velocity[2] = velocity_unit[2] * velocity_magnitude_heat;
    }

}



template<int dim>
void CellParticles<dim>::delete_useless_particles()
{
    sort(particle_delete_list_.begin(), particle_delete_list_.end());

    int i_particle_last = data_.size() - 1;
    for(int i_list = 0; i_list < particle_delete_list_.size(); i_list++)
    {
        //delete particles from end to begin (for the case: the last particle is to deleted)
        int i_particle = particle_delete_list_[particle_delete_list_.size() - 1 - i_list];
        if(i_particle == i_particle_last)
        {
            i_particle_last--;
        }
        else
        {
            //swap(data_[i_particle], data_[i_particle_last]);
            data_[i_particle] = data_[i_particle_last];
            i_particle_last--;
        }
    }
    data_.resize(data_.size() - particle_delete_list_.size());
    particle_delete_list_.clear();
}


/*
template<int dim>
void CellParticles<dim>::delete_useless_particles()
{
        for(int i_list = 0; i_list < particle_delete_list_.size(); i_list++)
        {
            //delete particles from end to begin (for the case: the last particle is to deleted)
            int i_particle = particle_delete_list_[particle_delete_list_.size() - 1 - i_list];
            data_.erase(data_.begin() + i_particle);
        }
        particle_delete_list_.clear();
}
*/

template class CellParticles<1>;
template class CellParticles<2>;
template class CellParticles<3>;
